Tomato Spotted Wilt Virus

In [1]:

    

using CSV, DelimitedFiles, Distances, Random, Pathogen, Plots
Random.seed!(5432)

# POPULATION INFOFORMATION
# Use CSV.jl for DataFrames I/O
# We know the types of the columns, so we'll manually specify those.
# * Individual IDs are `Int64`
# * X,Y coordinates are `Float64`s
risks = CSV.read(joinpath(@__DIR__, "plant_locations.csv"), types=[Int64; Float64; Float64])

# Will precalculate distances
distances = [euclidean([risks[i, :x]; risks[i, :y]], [risks[j, :x]; risks[j, :y]]) for i = 1:size(risks, 1), j = 1:size(risks, 1)]

pop = Population(risks, distances)


    

    
Out[1]:





Population object (n=520)

In [2]:

    

# OBSERVATIONS
# Use julia's included CSV interface for simple vector of observation times
raw_observations = readdlm(joinpath(@__DIR__, "infection_observations.csv"))[:]

# Create an `EventObservations` object with `Pathogen.jl`
obs = EventObservations{SI}(raw_observations)

# For performing inference we are going to set everything at or before t = 42 as being the starting state.
starting_states = [obs.infection[i] <= 42.0 ? State_I : State_S for i=1:obs.individuals]

# We will also set these observation times to -Inf
obs.infection[obs.infection .<= 42.0] .= -Inf

obs


    

    
Out[2]:





SI model observations (n=520)

In [3]:

    

# RISK FUNCTIONS
function _zero(params::Vector{Float64}, pop::Pathogen.Population, i::Int64)
  return 0.0
end

function _one(params::Vector{Float64}, pop::Pathogen.Population, i::Int64)
  return 1.0
end

function _powerlaw(params::Vector{Float64}, pop::Pathogen.Population, i::Int64, k::Int64)
  α = params[1]
  β = params[2]
  d = pop.distances[k, i]
  return α * (d^(-β))
end

rf = RiskFunctions{SI}(_zero, # sparks function - we will assume no exogenous transmissions and set this to zero
                       _one, # susceptibility function - we do not have individual level risk factor information to explore here, so will set to a constant 1
                       _powerlaw, # transmissability function - we will use a powerlaw (with intercept) kernel. This provides a spatial and non-spatial component to infection transmissions. This has 3 parameters.
                       _one) # infectivity function - we do not have individual level risk factor information to explore here, so will set to a constant 1


    

    
Out[3]:





SI model risk functions

In [4]:

    

rpriors = RiskPriors{SI}(UnivariateDistribution[], # empty `UnivariateDistribution` vector for all parameter-less functions
                         UnivariateDistribution[],
                         [Gamma(1.0, 0.5); Gamma(1.0, 1.0)], # Relatively uninformative priors with appropriate support
                         UnivariateDistribution[])

ee = EventExtents{SI}(14.0)

mcmc = MCMC(obs, ee, pop, starting_states, rf, rpriors)
start!(mcmc, attempts = 50000)
iterate!(mcmc, 50000, 2.0, event_batches = 20)


    

    




Initialization progress100%|████████████████████████████| Time: 0:09:32m17
MCMC progress100%|██████████████████████████████████████| Time: 17:50:5253m






    
Out[4]:





SI model MCMC with 1 chains

In [5]:

    

#using JLD2
#@save "mcmc.jld2" mcmc
#@load "mcmc.jld2" mcmc


    

In [20]:

    

gr(dpi=200)


    

    
Out[20]:





Plots.GRBackend()

In [21]:

    

p1 = plot(1:20:50001, mcmc.markov_chains[1].risk_parameters, yscale=:log10, size=(400,300))
png(p1, joinpath(@__DIR__, "trace.png"))


    

In [24]:

    

summary(mcmc, burnin=10000, thin=20)


    

    
Out[24]:




2 rows × 4 columns
parameter
mean
var
CI
String
Float64
Float64
Tuple…
1
κ₁
0.00265338
6.19118e-8
(0.00217468, 0.00313057)
2
κ₂
2.13986
0.0178767
(1.87603, 2.39948)

In [8]:

    

tnd = TNDistribution(mcmc, burnin=10000, thin=20)


    

    
Out[8]:





TransmissionNetworkDistribution Σexternal = 13.0 Σinternal = 314.0

In [27]:

    

p2 = plot(tnd, size=(400,300))
png(p2, joinpath(@__DIR__, "posterior_outdegree.png"))


    

In [28]:

    

p3 = plot(tnd, pop, size=(500,250));


    

In [29]:

    

png(p3, joinpath(@__DIR__, "posterior_tn.png"))


    

In [30]:

    

p4 = plot(mcmc.markov_chains[1].events[10000], State_S,
          linealpha=0.01, title="S", xguidefontsize=8, yguidefontsize=8,
          xtickfontsize=7, ytickfontsize=7, titlefontsize=11)
for i=10050:20:50000
  plot!(p4, mcmc.markov_chains[1].events[i], State_S, linealpha=0.02)
end

p5 = plot(mcmc.markov_chains[1].events[10000], State_I,
          linealpha=0.01, title="I", xguidefontsize=8, yguidefontsize=8, xtickfontsize=7, ytickfontsize=7, titlefontsize=11)
for i=10050:20:50000
  plot!(p5, mcmc.markov_chains[1].events[i], State_I, linealpha=0.02)
end
plot!(p5, obs, State_I, linecolor=:black, linewidth=1.5) # Show infection observations (day of prodrome)

l = @layout [a b]
combinedplots1 = plot(p4, p5, layout=l, link=:y, size=(800,400))
png(combinedplots1, joinpath(@__DIR__, "posterior_epi_curves.png"))